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These lecture notes are prepared as teaching material for an introductory course on Com- 
plex systems science for the use of master and PhD students in computer science and 
engineering. As such, they should be considered a work in progress and always under 
constant revision. 


The field of complex systems is wide and these notes cover only a portion of it, which 
the author believes to be relevant for engineering and computer science curricula. Most of 
this teaching material is taken from the documents cited in the bibliography, with revisions 
and amendments. However, the author of these notes is the sole responsible for errors and 
inaccuracies. 


“I am grateful to Prof. Silvana Bettelli Biolchini, who disclosed to me the beauty of science and intro- 
duced me to fractals, chaos and dynamical systems when I was a student at the high school. I would like 
to dedicate her these notes, hoping that she will enjoy reading them, remembering those pioneer times in 
which we had to wait a whole day to get the first fractal image on a computer monitor. 
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Introduction 


Complex system science (CSS) is the corpus of theories and methods that can help dealing 
with complex systems (CSs). Examples of complex systems are the brain, the society, the 
ecosystem, the cell, the ant colonies, the stock market. There is no formal definition of 
a complex system—and it is unlikely that a formal definition indeed exists—but we can 
provide a tentative, informal definition by saying that CSs are characterised by some of 
the following properties: 


- Composed of many elements 
- Nonlinear interactions 


- Non-trivial network topology (i.e. not a simple structure of relations such as a tree 
or a star) 


- Positive and negative feedbacks 

- Adaptive and evolvable 

- Robust 

- Levels of organisation (multiscaling) 


CSS has the objective of studying and modelling, and controlling CSs. It is interdisci- 
plinary and it involves many disciplines, such as mathematics, physics, computer science, 
biology, economy, philosophy, neurology and more. 


A peculiarity of complex systems is that local rules often propagate information in 
such a way that the whole system is subject to a dynamics that is not possible to under- 
stand solely on the basis of the system's constituents. Indeed, it is often said that some 
global behaviour emerge in a complex system from the interaction of its parts, governed 
by local rules. Despite its relevance for CSs, the concept of emergence remains somehow 
elusive. I quote below three definitions that provide a sufficiently clear framework for this 
phenomenon. 


e “The term emergence describes the onset of novel properties that arise when a higher 
level of complexity is formed from components of lower complexity, where these 
properties are not present.” (from P.L. Luisi, The emergence of life, Cambridge University 
Press, 2006) 


e “Emergence refers to the arising of novel and coherent structures, patterns, and 
properties during the process of self-organization in complex systems. Emergent 
phenomena are conceptualized as occurring on the macro level, in contrast to the 
micro-level components and processes out of which they arise.” (from J. Goldstein, 
Emergence as a Construct, Emergence 1 (1):49-72) 


e “Emergence refers to all the properties that we assign to a system that are really prop- 
erties of the relationship between a system and its environment.” (from Y. Bar-Yam, 
Concepts in Complex Systems, 2000. Available on the New England Complex Systems 
Institute website, http: //necsi.edu/guide/concepts/emergence.html) 


Strictly connected to emergence is self-organisation, which refers to phenomena in 
which some regular pattern emerges in a system, without a direct external control. Self- 
organization can be often observed in natural systems such as ant colonies. 

Besides emergence and self-organisation, another concept that pervades CSs is univer- 
sality, which refers to the fact that there exist behaviours appearing in many disparate 
systems: despite their differences at the microscopic level, some systems show the same 
behaviour at the macroscopic level. One of the goals of CSS is to identify the classes of 
local rules that generate such universal behaviours.! 


‘An very nice book illustrating “how complexity pervades biology”, and not only, is [25]. 
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1 Systems and models 


In the course of this path along some pillars of CSS, we will encounter terms that have to 
be defined. Often, the definition will be provided in an informal way. This is the case of 
the word system, which is usually referred to objects, things or entities (both concrete and 
abstract) for which it is possible to identify some kind of boundaries so that the system can 
be distinguished from its outside. (Of course, this is a generic definition) Therefore, once 
we have identified a system, it is possible to say what belongs to the system and what does 
not, i.e., what is outside of the system’s boundaries (its environment). The identification 
of a system depends upon the specific level of abstraction chosen and the viewpoint of 
the observer. Usually, we can also say that a system is composed of other entities, whose 
relationships are more strong than those with the outside of the system. 

The identification of systems is the first step into the study of the world around us; 
indeed, ancient philosophies and biology are full of discussions on systems and their clas- 
sification and taxonomies. However, another concept is needed to systems” science: the 
model of a system. A model is an abstract and schematic representation of a system. 
It is often a formal representation of the system. It has to be emphasised that a model 
represents a portion of the system, as it captures only some of its features. In addition, 
it requires an abstraction process, which involves simplification, aggregation and omission 
of details.? In general, we make a model of a system to (i) understand it and investigate 
some properties, (ii) control it, and (iii) make predictions on its future. Important classes 
of models are differential equations and difference equations. However, in general a model 
can be any system we decide to use to represent another system. For example, we may 
want to model a biological system (e.g., a cell) by means of a computer simulation or of a 
wet lab experiment. 

The study of general properties of systems and models is the main subject of systems 
theory. 


2 «Essentially, all models are wrong, but some are useful”, in Box, G.E.P., and Draper, N.R., (1987), 
Empirical Model Building and Response Surfaces, John Wiley & Sons, New York, NY, p. 424 


2 Dynamical systems 


A prominent class of systems concerns those systems that evolve in time. A model of a 
system in this class usually contains the rules that govern the time evolution of the state 
of the system. The system’s state at a given time is the set of values of the magnitudes 
that are relevant for the system and that are somehow measurable. 

It is commonly acknowledged to say that a system is dynamical if it can be described by 
models written in terms of differential or difference equations. Of course, dynamical sys- 
tems can be modelled by means of other kinds of models. An excellent book on dynamical 
systems is [26]. 

Some important concepts for dynamical systems are the following: 


e Trajectory: Sequence of states encountered during system’s evolution along time. 
The trajectory con end at a fixed point or at a cycle, or it can be more complex (see 
next sections). 


e Phase space: Space of the variables needed to characterise the system. 


e State space: Space of the states of the system. It coincides with the phase space 
when the set of variables needed to characterise the system coincides with the system’s 
state. 


Typical models for dynamical systems are ordinary differential equations and difference 
equations (a.k.a. iterative maps). 


2.1 Ordinary differential equations 


We assume that the reader is already familiar with differential equations and in this section 
we summarise the main points with respect to complex dynamical systems. 

Ordinary differential equations (ODE) provide a description of the evolution in time 
of the system’s state. For example, let us consider a naive model of population growth. 
Let N(t) be the population at time t of the species under study. The equation ruling its 
evolution may be: 


NN (1) 


where r > 0 is the growth rate. We can solve the equation analytically and obtain 
N (t) = Noe” (2) 


where No is the population at time t = 0. This equation tells that the population will 
grow exponentially. This is clearly not realistic, so a damping factor should be added 
to have a more realistic model. Indeed, when the population reaches a given level, some 
disadvantages usually arise due to resource competition (e.g., food and space). A cor- 
rection accounting for a linear dependence between N /N and N, leading to the logistic 


K/2 K 


equation:3 


; N 
N=rN|l- = 3 
r ( =) (3) 
where K is the carrying capacity. 


The equation can be studied analytically and the stability of its fixed points can be 
determined as illustrated in Box 1. Here, following [26], we study it graphically to have 
an idea of the qualitative behaviour of the system. We can plot N vs. N and study the 
equation as if it was a vector field. From Fig.1 we can observe that the equation has two 
fixed points, one at N = 0 and the other at N = K. These two fixed points are the only 
steady states of the system. The fixed point at 0 is unstable, while the one at K is stable. 
The qualitative behaviour of N(t) can also be derived from Fig.1 (see [26]). The study of 
the qualitative behaviour of Equation 3 tells us that, according to this model, a population 
with N(0) > 0 will grow fast at the beginning, then it will slow down its growing rate, until 
reaching a fixed point N = K from which it will not move, even in front of perturbations. 
Should we use this model to study a real system, we would also undertake a process of 
parameter tuning, in which the parameter K and r would be set to a value such that the 
behaviour fit some real data. Note that, since we have required to be r > 0, the behaviour 
of N(t) is qualitatively the same for every feasible value for r. 


ODE are widely applied and they can be way more complicated than the logistic equa- 
tion. Several branches of CSS make extensive use of differential equations, which are indeed 
a prominent class of models for complex systems. Another prominent class of models is 
that of difference equations, which we introduce in the next section. 


3The logistic equation is able to capture some simple growth phenomena, but it has several limitations 
and several variants have been proposed in the literature. 


Box 1: Stability of fixed points in differential equations 


An analytic way to determine the stability of a fixed point x* for a dynamical system defined 
by « = f(x) is the so-called linear stability analysis. The idea is to study the behaviour of a 
perturbation to x to see whether it decreases or increases (or keep the same value) in time. Let 
us consider a perturbation of the fixed point: n(t) = x(t) — x*. The behaviour in time of n(t) is 
given by n(t) = La — x*) = t. Thus 7 = f(x) = f(x* +n). Using a Taylor series expansion we 
have: 7) = f(a*) + nf’(a*) + O(n?). By observing that f(x*) = 0 and supposing that f’(2*) Æ 0, 
the quadratic term is negligible, so we have 1) = nf/(x“), which says that the perturbation size 
depends exponentially on t. Hence, the sign of f'(x*) tells us whether the perturbation is going 


to increase (f’(a*) > 0) or decrease (f’(x*) < 0). 


2.2 Difference equations 


Difference equations (DE) are used to model dynamical systems that update their state 
at discrete time steps. This class of models is widely used in the context of CS and it is 
particularly suitable to be studied in computer simulations. DE are of the form: 


x(t +1) = fiz); P, tl (4) 


where x, f and P can be vectors and represent the state of the system, the transition 
functions and the parameters of the model, respectively. When a system does not depend 
upon t it is named autonomous.’ For a detailed description of DE in complex systems 
modelling (and much more!), I suggest to read the beautiful book by Serra and Zanarini, 
Complex systems and cognitive processes [19]. 


Let us consider the case in which we need to model a system by means of a DE. Imagine 
we want to design a model for the evolution in time of the number of flies in a region.” 
The number of flies is measured ones a year for a number of years, so the data we have 
are inherently discrete in time. The main observation is that the number of flies in year 
t+ 1 depends on the number of eggs laid, which in turn depends upon the number of flies 
in year t. As a first attempt, let us try to model the dynamics of this system with a linear 
DE and suppose that for each fly in generation t there will be R flies in generation t + 1. 
In formulas: 


N(t+1)= RN; (5) 


The similarity between Eq. 5 and Eq.1 is apparent; indeed, they both model an expo- 
nential growth of the population. But there is more: Eq. 5 is the discrete-time version of 
Eq. 1. The qualitative behaviour of Eq. 5 can be easily studied by simulation. Once we 
have assigned a value to R and to N(0), we iterate the map for a number of steps and we 


4Note that this term is used also in other fields with a different meaning. 
5This example is taken from [12]. 


can trace the trajectory of the system in the state space. Obviously, this is not a proof, 
but it helps having an idea of the behaviour of the system. There are two ways to iterate 
a DE: the first one is called Cobweb method,” which is a graphical method; the second is 
an iterative numerical procedure, which can be conveniently implemented in a computer 
program. 

Let us assume we want to iterate numerically the equation. The key factor in this 
analysis is the value of R, which induces different kinds of behaviours. 


Case 1: Decay: When 0 < R < 1 it is easy to see that the population will decay in time. 


Case 2: Growth: When R > 1 the population increase exponentially (exponential 
growth). 


Case 3: Steady-state behaviour: When R = 1 the population stays at the same level. 


This simple analysis of the behaviour of the system depending on the values of the model 
parameters has the goal of emphasising two important points: 


1. once we have designed a model it is mandatory to study its behaviour as a function 
of its parameters; 


2. when we calibrate the parameters so as to fit real data we must pay attention to the 
regions in which the system will work—which may also vary in time! 


As in the previous subsection, we observe that Eq. 5 is a quite rough model of the real 
system. Therefore, we introduce a variant that makes the offspring per fly decrease as N (t) 
gets larger. The new equation is the so-called logistic map: 


N(t+1) = RN(t) - bN(t)? (6) 


where b5 0 rules the way in which the growth rate decreases with N(t). We can transform 


Eq. 6 by a change of variables: be x(t) = INO T For readability, we also change the 


R 
notation by writing t as a subscript. We then obtain: 


Te = Rx (1 — 1) (7) 


Here we assume that xy is the fraction of individuals in the population w.r.t. the maximal 
population size, so it holds x, € [0, 1]. 

This equation is a discrete version of Eq. 3 and we may expect that the system undergoes 
the same kinds of behaviour. As we will see, this is definitely not the case. To study the 
equation, we can draw z1 vs. x; (see Fig.2). As we can observe, the relation is quadratic, 
with the maximum equal to R/4 at x, = 0.5. We restrict the values of R to the range [0,4] 
so as to map z+ into itself. 


See [26, 12]. 
"We scale the number of flies by b/R. 


Box 2: Stability of fixed points in difference equations 


Given the iterative map 2,1 = f(x,) with a fixed point «*, we ask if the fixed point is stable 
and we apply the same idea as before. The perturbation is now Nn = £n — 1*; we have: x* + 
Ma = f(a* +n) = f(x*) + f'(2*)m + O(n?). Since 2* = f(x*) by definition, the perturbation 
changes in time according to the following difference equation: nn+1 = f (2%) + O(n). Thus, 
if f'(a%) 4 0, we have: mg © f’(x*)nn. Therefore, the perturbation increases exponentially in 
time if | f’(a*)| > 1 and decreases if | f’(a*)| < 1. 


By imposing x, = Rx,(1--x,) we find the fixed points of the map, which are xj = 0 and 
x3 = (R—1)/R (whenever R > 1). The stability of these fixed points can be studied as 
illustrated in Box 2. As expected, when R < 1 the only fixed point is xj, which is stable. 
For R > 1 both fixed points exist, of which the first is unstable and the second is stable. 
However, the situation is more complicated as it seems and it can be summarised as follows. 


Depending on the values of R we may observe different behaviours: 


e For small values of R, i.e. R < 1, the population always goes extinct (Fig.3). 
e For 1 < R < 3 the population grows and eventually reaches a fixed point (Fig.4). 


e For R > 3 the population oscillates (the attractor is a cycle). The period of this 
attractor depends on the value of R, starting from 2 and doubling while R increases 
(Fig.5 and Fig.6). 


Let us now consider this last case of oscillating behaviour. Period doubling occurs faster 
and faster: 


R Period 
Ri=3 2 
Ro = 3.449... 4 
Rs = 3.54409... 8 
R, = 3.5644... 16 
Rs = 3.568759... 32 
Ro = 3.569946... 00 


The value R, eventually converges to the value Ra, = 3.569946. The distance between two 
consecutive period doublings reduces as R increases: 


Seria BLOG (8) 
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Figure 2: Graphical representation of a numerical simulation of the function x,,, vs. £r, 
with R=2. 
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Figure 3: Graphical representation of a numerical simulation of Eq. 7 with R = 0.8 and 
xo E {0.25, 0.55, 0.85}. 
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Figure 4: Graphical representation of a numerical simulation of Eq. 7 with R = 2.2 and 
xo E (0.25, 0.55, 0.85}. 
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Figure 5: Graphical representation of a numerical simulation of Eq. 7 with R = 3.3 and 
Xy = 0.85. The period of the cyclic attractor is 2. 
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Figure 6: Graphical representation of a numerical simulation of Eq. 7 with R = 3.5 and 
x9 = 0.85. The period of the cyclic attractor is 4. 
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Figure 7: Graphical representation of a numerical simulation of Eq. 7 with R = 3.9 and 
Xo = 0.85. 
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Figure 8: Bifurcation diagram for Eq. 7. 


The limit is called the Feigenbaum's number. 


We may ask now what does happen for R > Ra. In fact, something literally strange 
happens: the attractor is no longer a cycle and it seems that the trajectory randomly 
moves across the state space (see an example in Fig.7). A random behaviour is excluded, 
since the model is deterministic. So what happens? We observe a chaotic behaviour, 
which will be the subject of the next section. 

Before moving to chaos, let us remark that the lesson learned from this example is that 
the behaviour of a system (more precisely, of the model of a system) may be dramatically 
different depending on the values of its parameters. A pictorial view of the different 
behaviours is provided by the bifurcation diagram, shown in Fig.8. 


2.3 Chaos 


The attractor of the system described by Eq. 7 for Ræ < R < 4 is neither a fixed point, nor 
a cycle, but it is an attractor of different nature which is called strange attractor. A strange 
attractor seems to be random, but in fact it is not. The trajectory of the system is confined 
in a specific area of the phase space, but it moves in a way that is almost unpredictable. 
The discovery of strange attractors is attributed to Lorenz, who was working on equations 
for weather forecast in 1963. In Fig.9 the Lorenz's attractor is shown. 
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Figure 9: Lorenz's attractor. 


The signature of chaos is the sensitive dependence on initial conditions: two nearby 
initial states produce trajectories that diverge exponentially in time. This property makes 
long term predictions very difficult, as doubling the precision in the forecast requires an 
increase of an order of magnitude in the precision of the measurements of the initial con- 
ditions. 

The sensitivity on initial conditions can be studied analytically (for example, by means 
of the Lyapunov exponent [19]) or it can be observed also by numerical simulations. This 
experiment on Eq. 7 is left to the reader as an exercise. 

Chaos is not random, so there should be ways to detect it, even in the cases in which 
we do not have access to the original model (or we are studying a real system). If we 
take a sample trajectory of a chaotic behaviour we can plot x444 vs. x, For a random 
map, such as 2,4; = Random(0,1) the plot would look like the one depicted in Fig.10. 
For the trajectory shown in Fig.7 we have instead the plot of Fig.11. Observe that Fig.11 
represents a quite regular relation, actually like the one depicted in Fig.1. A notable case 
of discovery of chaotic behaviour in a real system (heart) can be found in [9]. 


Since its first discovery, chaos has been observed in many natural and artificial systems, 
from weather to electronic circuits. It is important to observe that a chaos is not a synony- 
mous of complexity, as it will be discussed later in these notes. To know more on chaos, I 
suggest the book by Gleick [10]. 
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Figure 11: A chaotic map obtained from the trajectory in Fig.7. 
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Figure 12: The Cantor set (6 iterations of the generating procedure). 


Figure 13: Koch snowflake. Figure 14: The Julia set. 


2.4 Fractals 


Strange attractors are important, as they represent situations in which predictions in a 
system's dynamics are rather difficult. A notable property of strange attractors is that 
they can be geometrically characterised as fractals. Fractals are geometrical objects with 
the following properties:* 


e it has structure at arbitrarily small scales; 
e it is self-similar; 
e it has a non-integer dimension. 


Self-similarity in a fractal, means that the object is composed of smaller copies of itself. 
As an example of a fractal, let’s take the Cantor set. Given a segment, the Cantor set 
is obtained by recursively applying the following procedure: for each segment s in the 
set, divide s into three equal parts and delete the central one (Fig.12). The Cantor set is 
completely self-similar, because it is exclusively composed of smaller copies of itself. 

Other examples of fractals are the Koch curve (see Fig.13) and the Julia set (see Fig.14). 

By looking at these fractals, a paradox arises concerning their dimension: it seems not 
sound to define their dimension as the minimum number of coordinates needed to locate 
every point of the set, as we informally say for Euclidean objects such as segments, curves 
and planes. To account for a dimension of fractal objects, a more general definition is 


8See [16]. 
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needed. Here we introduce the similarity dimension [26], which we illustrate by means of 
some examples in the Euclidean geometry. Let's take a segment of length L: if we shrink 
it by 2 obtaining a segment of length Lı = L/2, then we'd need 2 segments of length Lı 
to obtain the first one of length L. If we shrink it by 3, then we need 3 smaller segments, 
and so on. Let’s take a square: if we shrink it by 2 in each direction, we need 4 smaller 
squares to obtain the first one; if we shrink the square by 3, then we need 9 other smaller 
squares, and so on. In summary, we have a scale factor r and m number of copies. The 
similarity dimension is defined as follows: 


logm 
-— (9) 
log r 
In the case of a segment, we have d lng = jog r 1. For a square: d pam E 2 
And for a cube? d = 3! Now, let's calculate the dimension of the Cantor set, which is 
recursively given by 2 copies of a reduction by 3: d = qam = oe = 0.63. The Koch curve 
depicted in Fig. 13 has dimension d = 224 = 1.26. 


log 3 
This notion of similarity dimension can be also extended to the cases in which we do 


not have m and r by construction. It is based on the idea of counting the boxes needed to 
cover a given object. 


Besides being interesting mathematical objects, fractals—and their extensions such as 
multi-fractals—provide also very useful models for real systems. Many growth phenomena 
are characterised by fractal properties, as well as urban development, market dynamics. 
In addition, fractals are applied into fields such as medical data analysis, graphics and 
material science. 
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3 Complexity: dynamical regimes between order and 
chaos 


In the previous sections we have seen that systems can behave either quite regularly, at- 
taining a fixed point or a cycling attractor, or they can be chaotic. However, our experience 
suggests us that another kind of behaviour exists, which is somehow between order and 
chaos. This is the realm of complexity, which characterises behaviours that are neither 
totally regular, nor completely chaotic. In a sentence, complexity lies at the edge between 
order and chaos [27]. CSS tries to characterise complexity, even to measure it. 

In this section we illustrate the principles of complexity by means of some prominent 
models. 


3.1 The Ising model 


The Ising model is a beautiful example of a system that can exhibit both ordered and 
chaotic behaviours and can illustrate some important properties of the region between 
order and chaos. The model is a classical example in statistical physics, but here we just 
consider its qualitative behaviour in a simplified version, which can anyway be subject to 
experimental validation through computer simulations.’ 

Let's take a lattice L x L of atoms characterised by a spin!?, which can be either up 
(+1) or down (—1). Each cell tends to align its spin according to the values of its first 4 
neighbours, so the ideal asymptotic states of the systems are composed either of all spins 
up or all spins down. 

The magnetisation of the system is defined as: 


1 
Kona (10) 


where s; is the spin of atom 7. 
The system tries to minimise its energy, defined as follows: 


E--N Jsisj (11) 
(b3) 


where J > 0 is a parameter accounting for the coupling between atoms and (i, 7) denotes 
the set of all neighbouring pairs—counted only once. 

The dynamics of the system is influenced by the temperature T: when T is low, the 
spins are likely to be flipped when the flip tends to align the atom with its neighbours; 
conversely, when T is high, spins tend to change even if their flip increases the energy. 


“The description of the Ising model we provide here is mostly taken from [24]. 
The chapter containing this subject can also be downloaded as a sample chapter at: 
http://press.princeton. edu/chapters/s9483. pdf. 

The original model aims at reproducing ferromagnetic phenomena in materials. 
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Figure 15: Ising lattice snapshot for T ~ Ty. 


The magnetisation is called the order parameter of the system and the temperature is the 
control parameter. 

The system can be studied analytically, but here we want to describe the main prop- 
erties of the system in a qualitative way and draw some conclusions from an experimental 
analysis based on Monte Carlo simulation, following Algorithm 1 which represents the 
so-called Metropolis algorithm. 


Algorithm 1 Monte Carlo simulation of a 2D Ising model. Adapted from [24]. 
while maximum number of iterations not reached do 
Choose a random atom s; 
Compute the energy change AF associated to the flip s; 4- —s; 
Generate a random number r in [0,1] with uniform distribution 
if r < e37 then 
Si — —Si 
end if 
end while 


The probability e *r is called the Boltzmann distribution and k is a constant, which 
we assume here, for simplicity, to be equal to 1. 

It is not difficult to simulate the system and observe the different dynamics of the 
magnetisation M (T) as a function of T. For low values of T, the steady state of the system 
will be composed of atoms mostly frozen at the same spin and (M(T)) will be close either 
to 1 or —1; for high values of T the spins will randomly flip and it will be (M(T)) = 0. 
The interesting part comes when we look at intermediate values of T, especially those 
close to a particular critical value T.. For values close to T.., the magnetisation fluctuates 
considerably around 0 and structures emerge in the lattice (see Fig. 15). At T = T, a phase 
transition occur: the system magnetisation undergoes a change in its possible steady state 
values, as depicted in Fig. 16. 

The main properties arising in general at the critical point are the following: 


"Ut can be proven that T, = 2.27. 
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Figure 16: Bifurcation diagram for the 2D Ising model obtained by a Monte Carlo simula- 
tion of a 30x30 lattice. Dashed vertical line at T = T, = 2.27 

- Large fluctuations 

- Long-range correlations 

- Percolation 

- Scale-free (1/f) 


These properties will be discussed in more detail in the following sections. Here we 
just remark that in the proximity of the critical point (of the control parameter), the 
order parameter has large fluctuations, as depicted in Fig. 17. This phenomenon suggests 
to use special care when we deal with system in a dynamical regime between order and 
chaos/randomness. 
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Figure 17: Variance of the average magnetisation values of Fig. 16. The shift w.r.t. Te 
is due to the approximation introduced by the numerical estimations (length of transient 
and of the time range for computing the averages). 


3.2 Phase transitions, symmetry breaking and behaviour changes 


In the previous subsection, we have informally introduced the fundamentals of phase tran- 
sitions by means of the Ising model. The phase transition observed is of the second order, 
i.e., the order parameter changes continuously as the control parameter is varied. Indeed, 
what changes is not the dynamics of the single instance of the spin system, but rather the 
possible future steady states of an ensemble of spin system instances. 

The analogous of this ensemble phenomenon for a single dynamical system is called a 
bifurcation, which denotes the situations in which a system can exhibit qualitative changes 
with respect to the variation of some control parameters. Examples are the propagation of 
fire, information or viruses depending on a propagation threshold, and the change in the 
structure of a system, such as the shift from vegetated to desert habitats. 

Mathematically, a bifurcation corresponds to the appearance—and disappearance—of 
new steady states, sometimes alternative. It is important to emphasise that in these lat- 
ter cases, the final attractor of the system is path dependent, i.e., it is sensitive to initial 
conditions and little fluctuations in the transient. A pictorial view of this phenomenon 
is provided in Figure 18, in which a ball can follow different paths depending on small 
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Figure 18: The epigenetic landscape. Picture taken from Slack, J.M.W., Conrad Hal 
Waddington: the last Renaissance biologist?, Nature, vol. 3, November 2002. 


fluctuations on its initial condition or along its trajectory, up to a point from which on it 
will be trapped in one of the possible future dynamics. 


Bifurcations can be illustrated by simple mathematical examples, such as the following 


differential equation: 
d 
The fixed points of the system are 0, +,/1, —y/1; of course depending on the sign of 
u we can have either one or three fixed points. The stability of a fixed point x* can be 
determined by checking the sign of In = —x*. Hence, for u < 0, the only fixed point 


is x* = 0, which is stable (see Figure 19). For u > O, the point 0 becomes unstable, whilst 
the new fixed points + are stable (see Figure 20). The properties of the steady states of 
the system can be represented by a bifurcation diagram, as depicted in Figure 21. 

The phenomenon of symmetry breaking can be understood also by means of a potential 
function ®,,(x), defined by the following equation: 


(x) 
dx 


Fila) == (13) 
If we suppose that the hypothesis of integrability of f,,(x) is valid —which is indeed our 

case—then we have: 
Dd, (1) = -Ër 4 gt (14) 


The function ©,,(x) is depicted for several values of u in Figure 22 (u = 1), Figure 23 
(u = 0.1) and Figure 22 (u = —1). The extrema of @,,(x) corresponds to the fixed points 
of f,(x) and the second derivative of @,,(x) characterises the stability of fixed points. 
(Note that the metaphor of a ball rolling down a smooth surface matches quite well the 
mathematical intuition) 
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Figure 19: f (x) for u = —1 
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Figure 21: Bifurcation diagram of the system described by Eq. 12. Picture taken from [24]. 


As we can observe, the system can undergo two kinds of steady-state behaviours, de- 
pending on the sign of u. Moreover, note that for u > 0 the system can reach either of two 
fixed points, so in a sense it has to decide between two alternatives. 

Symmetry breaking is also typical of complex systems exhibiting self-organisation— 
e.g., ants in foraging behaviour choosing between two paths of equal length. Therefore, we 
can understand the meaning of the critical point (in this case, u = 0): both in the case of a 
phase transition in systems composed of many elements and in mathematical bifurcations, 
the critical point separates two different kinds of system’s behaviour. It is also interesting 
to observe what happens for values of y close to the critical value: the potential is flat in 
the vicinity of 0, so a small perturbation might move the system far away from the initial 
point before it can return to the stable fixed point. This is a property typical of systems 
close to a critical phase: a perturbation takes very long to get dampened. This is also 
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Figure 22: @,,(x) for y = 1. 


Figure 25: Fire behaviour in forests at different density p. (a) p < pe; (b) p = pe; (c) 
p > pe. Picture taken from [24]. 


a consequence of the fact that in critical systems, information exchange takes part also 
between “distant” parts of the systems. 


3.3 Examples 


In general, we are interested in those situations in which the system can undergo qualitative 
changes as a function of some control parameters. In this section, we briefly review some 
notable examples. 


3.3.1 Forest-fire model 


A relevant property of many critical systems is percolation, i.e., the phenomenon in which 
information flows across the whole system. In situations in which percolation occurs, the 
system assumes a particular structure in which its elements form connected islands of 
activity. 

A case in point is the forest-fire model: the simplest model is provided by a lattice in 
which cells can be either occupied by a tree or empty. Cells are initially assigned by means 
of a Bernoulli distribution of parameter p. If all the trees of a given edge of the lattice are 


25 


O i 
#3 110 011 111 
000 
= 
001 ——=> 010 


(a) A BN with three nodes (b) State space 


Figure 26: An example of a BN with three nodes (a) and its corresponding state space 
under synchronous and deterministic update (b). The network has three attractors: two 
fixed points, (0,0,0) and (1,1, 1), and a cycle of period 2, {(0,0,1), (0,1,0)}. 


burned, then fire spreads across the lattice depending on the density of the trees. Below 
a critical value of p (pe = 0.59), the fire dies quickly, as trees are sparse; instead, above pe 
the fire reaches the whole area very fast because trees are dense. At p = p, there is a phase 
transition, in which the control parameter is p and the order parameter is the fraction of 
burned trees (see Figure 25). 

An important property of the system at p = p, is that the fire percolates across the 
system, i.e., it spans the whole system but leaves some “holes” at all scales; the set of 
burned tress at pe composes a fractal object. 

It is interesting to remark that at the critical point there is an abrupt transition in fire 
spreading: for p > pe the fire, i.e., information, pervades the whole system. This property 
is typical of many complex phenomena. 


3.3.2 Random Boolean networks 


Boolean networks (BNs) have been introduced by Kauffman [13, 14] as a genetic regula- 
tory network model. BNs have been proven to reproduce very important phenomena in 
genetics and they have also received considerable attention in the research communities 
on complex systems [2, 14]. A BN is a discrete-state and discrete-time dynamical system 
whose structure is defined by a directed graph of N nodes, each associated to a Boolean 
variable x;, i = 1,...,N, and a Boolean function f;(z;,,... Pia where K; is the num- 
ber of inputs of node i. The arguments of the Boolean function f; are the values of the 
nodes whose outgoing arcs are connected to node i (see Figure 26a). The state of the 
system at time t, t € N, is defined by the array of the N Boolean variable values at time 
t: s(t) = (x4(4),...,rn(t)). The most studied BN models are characterised by having 
a synchronous dynamics--i.e., nodes update their states at the same instant--and deter- 
ministic functions (see Figure 26b). However, many variants exist, including asynchronous 
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Figure 27: Critical line for RBNs, as a function of p and K. 


and probabilistic update rules [22]. 

BN models’ dynamics can be studied by means of usual dynamical systems methods, 
hence the usage of concepts such as state (or phase) space, trajectories, attractors and 
basins of attraction. BNs can exhibit complex dynamics and some special ensembles have 
been deeply investigated, such as that of Random BNs. Recent advances in this research 
field, along with efficient mathematical and experimental methods and tools for analysing 
BN dynamics, can be mainly found in works addressing issues in genetic regulatory net- 
works or investigating properties of BN models. 

A special category of BNs that has received particular attention is that of RBNs, which 
can capture relevant phenomena in genetic and cellular mechanisms and complex systems 
in general. RBNs are usually generated by choosing at random K inputs per node and 
by defining the Boolean functions by assigning to each entry of the truth tables a 1 with 
probability p and a O with probability 1 — p. Parameter p is called bias. Depending on the 
values of K and p the dynamics of RBNs is called either ordered or chaotic. In the first 
case, the majority of nodes in the attractor is frozen and any moderate-size perturbation is 
rapidly dampened and the network returns to its original attractor. Conversely, in chaotic 
dynamics, attractor cycles are very long and the system is extremely sensitive to small 
perturbations: slightly different initial states lead to divergent trajectories in the state 
space. 

RBNs temporal evolution undergo a second order phase transition between order and 
chaos, governed by the following relation between K and p: Ke = [2p.(1—p-.)]~', where the 
subscript c denotes the critical values [7] (see Figure 27). Networks along the critical line 
show equilibrium between robustness and adaptiveness |1]; for this reason they are supposed 
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Figure 28: A pictorial view of the sand pile model. Picture taken from [3]. 


to be plausible models of the living systems organisation. Recent results support the view 
that biological genetic regulatory networks operate close to the critical region [18, 23, 5]. 


3.4 Self-organised criticality 


Now a question arises as to why so many real systems are critical. The answer is not 
unique and still incomplete. However, we can say that some reasons for this phenomenon 
are connected with robustness and evolvability, others are to be found in the fact that 
criticality seems to be connected with a maximal information processing capability. A 
further reason is provided by the theory of self-organised criticality (SOC). 

SOC has been proposed by Bak [4, 3]. Succinctly, the theory states that many systems 
naturally evolve toward a critical attractor, i.e., an attractor at a phase transition. The 
interesting point here is that there is no need of tuning a control parameter, but rather 
it is the system that evolves toward this critical point. Self-organised critical systems are 
not in equilibrium and possess many degrees of freedom. Examples of SOC can be found 
in earthquakes, solar flares, evolutionary models, economic and social behaviours, just to 
mention some examples. 

SOC can be easily illustrated by the so-called sand pile model. Let us consider a device 
that periodically drops one sand grain at a time from the same position: grains fall down 
one over the other and create a pile (see Figure 28). While the pile grows, from time to time 
we observe avalanches, which are normally distributed and small on average. But, once a 
certain size of the pile is reached, avalanches are distributed according to a power law, i.e., 
there are avalanches of every intensity (see Figure 29). A system with this property is said 
to be scale-free. The pile can not grow over but it also can not stop the addition of grains, 
so it “reacts” by showing this behaviour which is a stable attractor. 
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Figure 29: Distribution of earthquakes of given energy. The probability P(s) of finding an 
earthquake of size s goes approximately as s~7. Picture taken from [3]. 


4 Information-theoretic measures in CSS 


In the study of complex systems we frequently resort to qualitative descriptions to identify 
some properties of a system. One reason is that some concepts, even though widely used 
since decades, are still elusive. Examples are the concepts of self-organisation and emer- 
gence, not to mention complexity itself. In this section, we will succinctly present some 
measures that can be usefully applied to quantify some properties of dynamical systems. 
The root of these measures is information theory.*? The link between complex systems and 
information theory is the fact that a computational process can be seen as the evolution 
in time of a dynamical system. Thus, some properties of the system can be conveniently 
studied by means of information-theoretic measures. 


4.1 Information entropy 


The first notion to be considered is (information) entropy, introduced by Shannon in 1948.73 

Let us consider a simple system of which we observe the state at a given time. The 
observation can be modelled as a random variable X which can assume values from a 
finite and discrete domain X. If the observation is x € Y, which has a probability P(x), 
then the information content of the observation is measured as EFE = —log2P(x).* An 
improbable observation conveys more information than one associated to high probability. 


12For a detailed and profound discussion on these subjects, I suggest to read the Information theory for 
complex systems lecture notes by Prof. Kristian Lindgren, http: //studycas.com/c/courses/it, visited 
on April 2014. For a quick overview on the subject, see [17]. 

13 The measure is also called the Shannon entropy. 

14The usual unity measure in information theory is the bit, so logarithms are in base 2. Hereafter, we 
will omit the base in the mathematical notation. 
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Figure 30: Entropy values for a Bernoulli distributed stochastic variable. 


We can characterise the system by averaging over its possible outcomes: 


H(X) =- Y P(e) logP(x) (15) 


TEX 


In the definition of A(X) we assume 0 log0 = 0. 


Intuitively, H(X) measures the degree of randomness of the process. For example, let's 
consider the case of a pure random binary variable: we have ¥ = {0,1} and P(0) = P(1) = 
0.5. Hence, H(X) = —2(0.5 log(0.5)) = 1. In general, for n symbols appearing with equal 
probability, we have H(X) = log(n), which is the maximum value. Conversely, if P(0) ~ 1 
(so P(1) = 1— P(0) = 0), the entropy is quite small, as it is very likely that X will assume 
value 0 (see Figure 30). Therefore, a high entropy characterises systems that show string 
tendency to disorder, whilst it is low for ordered systems. 


Example: entropy of the Ising model 


Let's consider the Ising model, as defined in Section 3.1. An estimation of its information 
content can be given by computing the average atom entropy. For atom S;, the probability 
of occurring with spin —1, is estimated by computing the frequency of occurrence of spin 
—1 in the series of states collected. As we can observe in Figure 31, the average entropy is 
low for low temperature values and high for high temperatures. Note that the relation is 
nonlinear, with a sharp increase at the phase transition. 
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Figure 31: Average entropy for a 30 x 30 Ising model, as a function of T. 


4.2 Related information measures 


Entropy is the base for computing other important information-theoretic measures. 

First of all, this definition can be extended to blocks of symbols of length n in s, so as 
to take into account also correlations among symbols. This leads to the definition of the 
block entropy of length n: 

H,-- ` P(sn) log P(sn) 
sn ES 
The entropy excess’ is defined as the difference between block entropies of length n and 
n — 1 and estimates the information required to predict the n-th symbol conditioned to 
the observation of n — 1 preceding symbols. In formulas: 


hn = AH, = H, ma HA 


We can extend this process to the second derivative (in discrete domains) and obtain the 
correlation information from length n: 


kn = A?H(n) = —H(n)+2H(n—1) — H(n=2),n>3 


Intuitively, the peaks of k, identify significant block regularities, i.e. maximum gain in 
information for specific block lengths. 


15Not to be confused with the excess entropy [21], which is defined for n > oo. 
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Also the mutual information I(X, Y) between random variables X and Y is defined in 
terms of entropies and estimates the average information one gains about Y after the 
observation of X, and viceversa: 


I(X,Y) = H(X) + H(Y) - H(X,Y) 


where H(X, Y) denotes the conjunct entropy of X and Y. 
A generalisation of the mutual information is the so-called integration (a.k.a. intrinsic 
information or multi-information): 


I(X1: X25 ..3X4) = Y A(X) — 1(X1, Xo, ..., Xp) (16) 


4.3 Algorithmic complexity 


The notion of entropy suggests to focus on the amount of disorder of a system and it can be 
computed by estimating the probabilities of occurrence of the symbols produced. A related 
measure does not require to know this probability distribution and relies on computability 
principles. The algorithmic complexity of a sequence of symbols refers to the length in 
bits of the minimal program that generates the sequence, without input. This measure 
is known as Kolmogorov complexity (KC). For example, the sequence 0101010101...01 has 
low KC, while a random generated sequence of bits has the highest value of KC. 

The KC can not be computed exactly, but it can be estimated, for example by means 
of compression algorithms. Note that a high compression ratio means that the sequence 
contains high regularity, while low compression ratio might be the signal of randomness. 
Of course, these estimations should be taken cum grano salis, because they depend on the 
compression algorithm used. 


4.4 Statistical complexity 


We could be satisfied with the previous measures to study the system, but in fact they do 
not capture the intuitive notion of complexity, as a property characterising systems at the 
edge of order and disorder. A random sequence of symbols 0 and 1 should be evaluated 
with low complexity, rather than its maximum. To overcome this problem, several mea- 
sures as been proposed [17] to account for statistical complexity (SC), i.e., the algorithmic 
complexity of a program that reproduces the statistical properties of a system. In this 
light, the SC of a random sequence is low. 


Among various measures of SC, we mention a simple yet effective one, which is called 
LMC complexity, by the name of its inventors [15]. The idea is rather simple: if we want 
the SC of a system to be high in intermediate regions between order and disorder, we 
can define a measure as the product of a measure that increases with disorder and another 
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Figure 32: Average disequilibrium for a 30 x 30 Ising model, as a function of T. 


which decreases with it. We know that entropy measures disorder, this is the former factor. 
The second can be the disequilibrium: 


D(x) = Y (Pto) - aa) (17) 


TEX 


The disequilibrium estimates the extent to which a system exhibits patterns far from 
equidistribution. For example, if the trajectory of a system is composed of only few of the 
possible states (e.g., a short cyclic attractor), then it has a high disequilibrium. 

The LMC complexity is defined as: LMC(X) = H(X)- D(X). 


Example: LMC complexity of the Ising model 


Let us computed the LMC complexity for the Ising model, as previously done for the 
entropy. Results for the disequilibrium and LMC complexity are depicted in Figure 32 
and Figure 33, respectively. As we can observe, the shape of the LMC complexity shows a 
remarkable peak at the phase transition, exactly where criticality resides. 
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Figure 33: Average LMC complexity for a 30 x 30 Ising model, as a function of T. 


Further notes on complexity measures 


Measuring complexity is still an elusive task and it is tightly entangled with the formal 
characterisation of emergence and self-organisation. The discussion of this extremely rele- 
vant topic is not in the scope of these lectures, so I just suggest some references to know 
more on the subject (8, 11, 6, 20, 17]. 
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